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We study the 6-dimensional dynamics - position and orientation - of a large sphere advected by a 
turbulent flow. The movement of the sphere is recorded with 2 high-speed cameras. Its orientation 
is tracked using a novel, efficient algorithm; it is based on the identification of possible orientation 
'candidates' at each time step, with the dynamics later obtained from maximization of a likelihood 
function. Analysis of the resulting linear and angular velocities and accelerations reveal a surprising 
intermittency for an object whose size lies in the integral range, close to the integral scale of the 
underlying turbulent flow. 



I. INTRODUCTION 

The advent of resolved Lagrangian measurements has 
helped understand the dynamics of turbulence from the 
point of view of fluid particles [ ]. In the experiments, 
solid tracers are followed in lieu of fluid particles, which 
naturally raises the question of the understanding of the 
dynamics of finite size objects in turbulent flows. It is 
a subclass of the issue of the dynamics of inertial parti- 
cles, i.e. particles who have inertia with respect to the 
fluid motions, either because their density differs from 
that of the fluid or because their spatial extent cannot 
be ignored. If the particles are quite small compared 
to the smallest fluid motion (the Kolmogorov dissipative 
length scale 77), arguments show that they behave as trac- 
ers of fluid motions. Observations have revealed a very 
intense intermittency in the motion of fluid tracers [2, 3]. 
They experience very strong accelerations, with a prob- 
ability distribution which displays stretched exponential 
tails [4]. 

When the diameter D of the advected particles is of the 
order of, or larger than 77, their equation of motion is not 
known (see, however [5-7]). We restrict our discussion to 
neutrally buoyant spheres. Several recent studies [ ] 
have shown that the acceleration statistics of such inertial 
particles does not gently reduce to a Gaussian behavior 
as their diameter increases. It is an important feature 
because the characterization of forces acting on an object 
advected by a turbulent flow has many applications in 
engineering (from mixing issues in industrial processes 
to dispersion in the oceans or in the atmosphere). 

The study reported here takes a leap forward in size 
and considers the motion of a neutrally buoyant sphere 
with diameter D of the order of the integral scale L- int (the 
scale at which energy is fed into the flow). In addition, we 
aim at resolving the six degrees of freedom of the particle 
dynamics, i.e. the goal is to obtain a tracking in time and 
space of the particle's linear and angular motions. This 
allows the study of the forces and torques acting on a 
(large) inertial particle. 



The tracking of the particle position in space can be 
carried out by using methods already developed and suc- 
cessfully tested for small particles [ ]. In comparison, 
following the orientation of the particle is much more 
challenging, both because of the specifics of angular vari- 
ables, and of specific algorithmic requirements. Previ- 
ous studies on tracking the orientation have focused on 
measuring one component of the angular velocity. The 
particles [13] used for this purpose are transparent, and 
contain an embedded mirror and a diameter of less than 
50/im, which is of the order of the Kolmogorov length 
scale, 77. In the experiments reported in [13], the angular 
velocity at a point neither the translation nor the angular 
motion could be tracked for very long. The principle used 
here is completely different: it consists simply in painting 
the particle with a suitable layout, and in retrieving its 
orientation. For algorithmic efficiency (and robustness) 
this is not done step by step but for the entire trajectory 
using a global path extraction. 

The text below is organized as follows: we first present 
the experimental setup and recall some important fea- 
tures of the orientation algebra in 3D. We then describe 
how the particle images are extracted from the movie 
images, and compared to synthetic images with arbi- 
trary orientations. Possible candidates are identified and 
then assembled into an orientation time series using a 
global maximization of a likelihood function. Finally, we 
present some results concerning the particle dynamics. 



II. BASICS 
A. Experimental Setup 

A turbulent flow is generated in the gap between 2 
counter-rotating impellers of radius R = 10 cm fit- 
ted with straight blades 1 cm in height. The flow do- 
main in between the impeller has characteristic lengths 
H = 2R = 20 cm and the working fluid is a water-glycerol 
mixture, whose density can be finely tuned. In order to 
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FIG. 1: Sketch of the experimental setup: a) image of the 
von Karman mixer; b) sketch of the camera arrangement; c) 
textured sphere for different orientations. 



be able to perform direct optical measurements, the con- 
tainer is build with flat Plexiglas (Poly [methyl methacry- 
late]) side walls, so that the cross section of the vessel 
is square. This type of von Karman swirling flow has 
been used extensively in the past for the study of fully 
developed turbulence [ ] ; its local characteristics approx- 
imate homogeneous turbulence in its center, although it 
is known to have a large scale anisotropy [14, 15]. A 
sketch of the setup is provided in Fig. 1 - further details 
about the flow turbulence are given later in section IV. 

A white, Poly Amid sphere with diameter D = 18 mm 
(accuracy 0.01 mm, Marteau & Lemarie, France) moves 
and rotates in the turbulent flow. It is neutrally buoyant 
in the fluid - whose density is adjusted to that of the 
particle p p = 1.14 g.cm -3 by addition of glycerol to wa- 
ter. The density mismatch, measured from sedimentation 
speeds, is found to be less than Ap / p = 10 -4 . The parti- 
cle is textured black and white by hand using either black 
nail polish or a black-ink permanent marker. Its mo- 
tion is tracked using 2 high-speed video cameras (Phan- 
tom VI 2, Vision Research) which record synchronously 
2 views at approximately 90 degree. The flow is illumi- 
nated by high power LEDs and sequences of 8 bit gray 
scale images are recorded at a rate of 600 frames per 
second. 

Both cameras observe the measurement region with 
a resolution of 650 x 650 pixels, covering a volume of 
15 x 15 x 15 [cm 3 ]. Hence, the particle diameter is 70 — 90 
pixels. In the choice of the particle texture, several fea- 
tures have to be considered: 

- a single view should correspond to a unique orientation. 

- illumination inhomogeneities may cause regions to look 
similar in the camera images. Optically resembling views 
should correspond to clearly distinct orientations. 

- the cameras are grayscale so the texture has to be black 



and white. 

- the number of black and white pixel should be approx- 
imately the same in every possible view. 

In our configuration, the camera can store on the or- 
der of 15,000 frames in on-board memory, thus limiting 
the duration of continuous tracks. The movies are down- 
loaded to a PC, waiting to be processed. The processing 
is done on a gaming PC with a state of the art graphics 
card. Algorithm development and code test is done on 
an Apple Macbook Pro. The code is written in Matlab 
2009a using the image and signal processing toolboxes as 
well as the Psychtoolbox extension [16, 17] which provide 
OpenGL wrappers for Matlab. 



B. Angular Variables 

The parametrization of an angular position in 3D space 
causes a number of difficulties which are briefly addressed 
in this section (see e.g. [18-20] for a more complete pre- 
sentation). One of them is caused by the degeneracy of 
the axes of rotation for certain orientations (the 'gim- 
bal lock' problem). Another is the choice of a suitable 
measure of distance between two orientations. 



1. Describing Orientations 

As stated by the Euler rotation theorem, 3 parameters 
are needed to describe any rotation in 3D. We use here 
Euler angles with the Tait-Bryan convention as shown 
in Fig. 2. In the transformation from Lab to Particle 
coordinate system (CS), we first apply a rotation around 
the z— axis of angle 9 Z1 followed by a rotation around the 
intermediate y— axis of angle 6 y and last a rotation of 
angle 6 X around the new x— axis. The rotations work on 
the object using a right handed coordinate system and 
right handed direction of rotation. We will denote an 
orientation triplet by an underscore, e.g. 6, in order to 
distinguish them from vectors (which are typeset in bold 
font, e.g. uS). 

Lab Particle 

W W R(e*) 

FIG. 2: Tait-Bryan rotation sequence describing the sphere's 
orientation. 

The orientation of the object is fully described by an 
orthogonal 3x3 matrix R, obtained from the composition 
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of the 3 elementary rotations: 

C9y C0 Z —COy S6y 

s0 x sO y c0 z + c0 x s0 z —s0 x s6 y s6 z + cO x c6 z —s6 x c6 y 
—cO x sO y cO z + sO x sO z cO x sO y sO z + s6 x cO z cO x cO y 

a)" 

with c- = cos (•) and s- = sin (•). Consequently, from any 
rotation matrix the 3 Euler angles can be extracted using 



9 — {@x, 0y,0z} 



'atan2(-E 12 ,E n ) N 

asin(E 13 ) 
V atan2(-E 23 ,a 3 3)y 



(2) 



enforcing X ,0 Z G [0, 2tt[ and 9 y G [— 7r/2, tt/2]. However 
this choice is not unique because there is a second triplet 
with H (9 X + 7r, sign(6g - 7T — y ,0 z + n) =J^(0 x ,9 y ,0 z ). 
Needless to say, multiples of 2ir can be added to each 
angle. An important practical consequence is that even 
for small changes in orientation the difference between 2 
Euler angle triplets, 0_ x and nas formally 4 possible 
results. 

The curvilinear coordinate is related to the angular 
velocity, cj p , (in the particle frame) by 



1 o 

c6 x 
s9 r c9 r c9, 



s0 y 

s6t. c8„ 



dt 



(3) 



H(^, 0y) 



dt 



For cos(0 y ) ~ 0, the determinant of the matrix H, 
det(H), vanishes and its inverse is not defined. In other 
words, finite body rotations need infinite change in the 
Euler angles. This singularity is called a gimbal lock 
and is a well-known problem in robotics and aerospace 
engineering. Geometrically, the second rotation turns 
the first axis parallel to the third axis of rotation, and 
the rotation loses 2 degrees of freedom. Unfortunately 
gimbal locks cannot be avoided by a wise choice of 
representation. 

One then needs to define a distance between 2 arbi- 
trary orientations, immune to this type of singularity. A 
natural distance between two arbitrary orientation ma- 
trixes, A and B, is 



(at) 



(4) 



Tr \JA - i|) (A - S) J = 6 - 2 Tr 

= 4(1 - cos (0)) 

using the A A T = B B T = 1 and that AB T is a rotation 
matrix with the eigenvalues l,e^,e - ^. The distance 
is thus a growing function of <p. We measure here the 
distance between two rotation matrices by: 



d (A, B) = acos 
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(5) 



Because it works directly on the orientation matrices it 
is neither sensitive to gimbal locks nor to the choice of 
the representation and thus an important tool in our 
algorithm. It should be noted that d(A,B) is the angle 
of the rotation which turned the orientation from A to B. 

In the search of the particle orientation, one last in- 
convenience of Euler angles is that they are not locally 
orthogonal, in the sense that 



d({0 x , By, 6 Z }, {9 X + A9 X , 6 y + A0 y , 6 Z + A^}) 2 a 
A8 2 x + AOl + A0 2 z + 2A9 X ■ A9 Z ■ em(0 v ) 



(6) 



for A small. As a consequence, a uniform spacing of the 
Euler angles in x: y ,0 z does not sample the space of 
possible orientations in an optimal way. In particular 
near gimbal locks, the sampling rate would be higher at 
no higher accuracy. The so-called Lattman angles [ ] 



9^0-} = {0 X + 9 Z , 0y, X — Z } 

fulfill local orthogonality since they verify 

d{{0 + , 0, #_}, {0 + + A# + , + A0, 0- + A#_}) 2 « 
A0 2 + (1 + sin 0) /2 + A0 2 + A0 2 _ (1 - sin 0) /2 . 



(7) 



(8) 



As they are locally orthogonal, it is sufficient 
for sampling purposes to keep A0^_(1 + sin#)/2, 
A/9 2 and A0 2 _(1 - sin0)/2 constant. After 
a constant sampling of N values of with 
^Latt = A0 = tt/(N — 1), the stepping in 0+ and 
6L can be computed with A0+(0) = AL a tt/ sm (f + f ) 
and A6L((9) = A Latt /sin (| - f). It should be empha- 
sized, that 0- G [0,27r[ whereas 0+ G [0,47r[. Lattman 
angles enable us to sample the set of orientations in an 
optimal way. 

Finally, in several instances it is convenient to describe 
a rotation by the direction of an axis n about which the 
systems is rotated by an amount <p. The corresponding 
rotation matrix can be computed using the Rodrigues 
Formula [18, 20] 

S(n,0) = 

c<j) + n x A —n z s(j) + n x n y A n y s(j) + n x n z A 
n z scj) + n x n y A ccj) + n 2 A —n x s(j)-\-nyn z A 



-n y s(j) + n x n z A n x sq 



■ riojUzA 



■ niA 



with A = (1 



(9) 



Eq. (9) also allows us to extract the axis, h , and the 
angle, <j) from any arbitrary rotation matrix. As a 
result, changing the coordinate system or changing the 
representation of rotation can be done by expressing the 
orientation in its matrix form, applying the transforma- 
tion which changes the CS and extracting the desired 
representation. 
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2. Angular Velocity and Acceleration 

Angular velocity and acceleration are often obtained 
by direct differentiation of a time-series of Euler angles, 
e.g. using Eq. (3). However, it is possible to obtain the 
angular velocity in the particle frame directly from the 
matrices. This technique is not sensitive to Gimbal locks 
because of the uniqueness of the orientation matrices. 

Let G^\y iZ be the particle CS at time step fc, whereas 
the fixed lab CS is . For two time-steps, k and 

k + m, we know the corresponding orientation matrices 
which rotate the particle: 



a(#fc+ m ) 

T 



-'x,y,z T 

L g(0fc + m )^ 

°x,y,z ^ 



^x,y,z 
D P,fc+ra 

D P,fc+ra 
"x,y,z 



in which the matrix T is the change in orientation, in 
other words the matrix representation of the discrete an- 
gular velocity (for a given time difference). The change 
is with respect to the particle CS at time k: T ex- 

pressed in the axis- angle convention (see Eq. (9)) returns 
a direction vector, n , of length unity and an angle, Acj) 
(meaning that between times k and k + m the particles 
has rotated an angle around the vector n ). The 
time difference, At, between the steps is a function of m. 
Therefore an estimator of angular velocity is 



J*(t(k)) = 



Acf) 
At(m) 



(10) 



Averaging over several separations, m, returns the 
angular velocity in the particle frame without a prior 
unwrapping nor problems near gimbal locks. The angular 
velocity with respect to the lab CS is defined as 



(11) 



The angular acceleration in either particle or lab frame 
is defined as 



OL 



L/P 



,L/P 



(12) 



In practice, it is obtained from a convolution of the angu- 
lar velocity time series with the derivative of a gaussian 
kernel. This technique has proved to be efficient in re- 
moving noise [ ]. 



III. TRACKING 

A. Position 

Although the identification of a large sphere from the 
camera images causes no particular conceptual difficulty, 



the fact that the sphere is textured raises some practi- 
cal issues. A simple thresholding returns only either the 
white or the black part of the particle. Reflections from 
the impellers continuously change the background, and 
small impurities in the flow and possible bubbles add 
sharp gradient noise to the images. Furthermore, the il- 
lumination of the flow is not perfectly uniform, and thus, 
shadows as well as reflections occur. 

For each movie and for each camera, we compute the 
background view as the average of an equally distributed 
subset of its images. For each frame we then subtract 
the background and perform a Difference of Gaussians 
blob detection. The threshold is adjusted by hand for 
each camera and light arrangement. Mat lab's Image- 
processing toolbox is used to identify blobs with a round 
shape and a diameter close to that of the particle. Shad- 
ows, bubbles, and reflections might be found during blob 
detection because of their sharp separation from the 
background, but they are of uniform texture and hence 
characterized by a small value of the variance of light in- 
tensity across the blob. The blob with highest variance 
and closest resemblance to a sphere is considered to be 
the particle. The precise position of the particle is refined 
using a Canny edge detection in a tight region around the 
blob. For each time step we record the position, (x,y), 
of the particle on the image in pixels plus its diameter, 
2 r, and the deviation from the spherical shape as an er- 
ror estimator. Since only one particle is placed into the 
flow, the track assembly is straight forward. The algo- 
rithm may temporarily loose the particle for short times 
(because of bad light reflection, blurs, ...); this is compen- 
sated by the large oversampling and gaps of less than 5 
frames are interpolated to obtain longer tracks. Outliers 
are be identified using a least square spline and replaced 
by an interpolation. 

Tsai's camera model and calibration technique [22] is 
used to project the 2D positions into 3D. The calibration 
of the cameras contains the position of the camera plus 
its rotation with respect to the lab CS, which is needed 
later for the orientation processing. 



B. Orientation 



The algorithm used to process the camera images and 
obtain a time series of orientations (and angular veloc- 
ities) can be split into 3 parts: (i) by comparison of 
the sphere's picture with synthetic images, the algorithm 
identifies a set of possible orientations; (ii) from the set 
of possible candidates at successive instants, a Flow algo- 
rithm identifies a likely time series; (iii) a post-treatment 
adjusts remaining ambiguities. These steps are described 
in details in this section. 
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1. Candidate Finding 

a. Synthetic images. A first step is to obtain a 2D 
projection, S(0), of a sphere with known texture and 
size at an arbitrary orientation, 0. This rendering is 
achieved using OpenGL, via the Psychtoolbox extensions 
for Matlab - for a disk image of about 60 pixels, the 
algorithm can render several thousand orientations per 
second (see Fig. 3 for an illustration). 

b. Texture extraction. Once the particle position 
and diameter are known, one extracts a disk subset of 
the image, centered on the particle, G. In a first step 
the contrast is adjusted such that the global histogram 
of intensity contains at least b percent of black and w 
percent of white pixels (the algorithm only takes into 
account the disk / particle region in G). The adjustable 
parameters 6, w are fixed to b = w ~ 30% which is the 
minimum amount of black/ white pixel in an arbitrary 
orientation. In a second step, the image is thresholded 
using Otsu's method [23] for the global histogram as well 
as for 2 moving regions. The thresholded image, I, is 
adjusted such that pixels outside the particle / disk are 
set to whereas black is —1 and white +1. These steps 
are shown in Fig. 4. 

c. Comparison, possible orientations. The image I 
(with diameter 2r) obtained as above is ready for com- 
parison with synthetic images. The resemblance to a 
rendered image S(0) with orientation is estimated us- 
ing the projection 

T ^ = l + ^EEky^-)> (is) 

i 3 




FIG. 3: Synthetic 2D projections of the particle for a range of 
orientations, using OpenGL. A camera image of the moving 
particle is shown in the upper left corner (contrast enhanced; 
note the driving disks on either side). 



which is ratio of the number of correct pixels to the total 
number of pixels. 

At this point we note that the computational cost of 
directly comparing an image I to synthetic ones S(0) cov- 
ering the set of possible orientation {6} scales roughly as 
^Latf wnere ^Latt is the grid spacing in the orientation 
space. There is also the additional difficulty that the 
particle apparent diameter changes slightly as the sphere 
moves in the flows. For efficiency and physical correct- 
ness, we use the following strategy: instead of finding at 
any time step the best images, we identify a set of possi- 
ble candidates for all time steps and then extract globally 
the time series of orientations. 

First we render images, S({# coarse }), covering all pos- 
sible orientations with a coarse grid - in practice ALatt ~ 
12°. Lattman angles are locally orthogonal and thus 
more efficient in creating such grids. The size of the 
rendered images is fixed to approximately one half of the 
particle real diameter. Since their size does not change, 
these images are kept in the computer memory and do 
not need to be recomputed for every new image. 

The thresholded particle image, I, is then resized to 
the size of the renderings, I , and compared to all 

=coarse 

synthetic images, S({# coarse }) as shown in Fig. 3 using 
Eq. (13). All angles with > max(T(I coarse , R oarse })) - 
^coarse are considered to be possible orientations (PO). 
Here S coarse is an arbitrary thresholding value, with in- 
spection showing that a value equal to 0.1 gives good 
results. 

Experience shows that the identified POs usually cover 
several broad classes. They are thus separated into 
groups of images whose orientations differ by less than 
a rough threshold, approximately 30 — 45°. For each 
group, synthetic images are further added using a fine 
grid spacing, Afi ne = 3° (at this point 'bad' images may 
cause the code to runaway; they are dropped and the 
code advances to the next time step). The PO images 
are then rendered in real size and compared (using the 
projection T) to the image I. For each group, the code 
returns the final best guess, i.e. the orientation with the 
maximum resamblance, thus drawing a list of candidates, 
see Fig. 5 for an example of a particle with its correspond- 
ing candidates. 

2. Track Assembly 

After identifying the candidates for each time step, the 
most likely orientation for each time step has to be deter- 
mined. However, the candidate with the highest count 
of correct pixels is not necessarily the best choice. Al- 
though counterintuitive, the direct use of 2 cameras see- 
ing the particle at different angles does not simplify the 
problem, because in the case of a bad image, one cam- 
era falsifies the choice of the candidates found by the 
other camera. Moreover, gimbal locks prevent the use 
of a predictor-corrector scheme for the prediction of the 
orientation. However, the norm of angular velocity is as- 
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FIG. 4: Texture extraction and comparison with a synthetic image. The resemblance between the image I and the synthetic 
projection S at angle is estimated using Eq. (13). 
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FIG. 5: Particle camera image (left) and corresponding can- 
didates, after analysis of the possible orientations (steps la-c 
described in the text). 
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FIG. 6: Sketch of a graph connecting the possible candidates 
using the cost function C (cf. Eq. (15)). In this example, no 
candidate could be identified at time setp 2. 



sumed to be smooth and we search the time series which 
globally minimizes the sum s ^2 t £(£) along the time series 
of the so called direct neighbor distance function: 

£{t) = \<jj(0(t),0{t + At)) \ =d(0{t),0(t + At)) /At . 

(14) 

A direct neighbor is the next valid time step at t + At. 
The distance between 2 orientations does not depend on 
the representation, ensuring the robustness of the algo- 
rithm even at gimbal locks. Minimizing £(t) is only 
meaningful for small changes in orientation between two 
time steps, another requirement for high (over) sampling 
rates. 

Flow algorithms are highly efficient in finding a global 
optimum for a discrete set of candidates. The following 
is done for each camera without considering the extra 
information from the second camera. In a first step we 
remove all candidates with a resemblance T < 5 qua ii ty - 
in practice <s qU aiity = 0.5. Then a directed graph is built 
which connects all candidates at time step t with all their 
direct neighbors at the non-empty time step t + At. The 
cost function is chosen such that it takes into account 
both the change in orientation and the quality of the 
matching: 

C({6 A ,T A },{6 B ,T B })=d{9 A ,9 B ) 2 ~ T ^~ TB , 

(15) 

with {9 a ,Ta} a candidate at time t and {0 b ,Tb} a di- 
rectly neighboring candidate at t + At. 



A Dijkstra path finding algorithm returns the sequence 
of candidates having a global minimum of the total 
cost, i.e. the global minimum of change of orientation 
(weighted by the image quality) (cf. Fig. 5). In most 
cases this algorithm returns directly the time series of 
absolute orientation. Nevertheless, bad images introduce 
false candidates forcing the path finding algorithm to 
take a different, non-physical path. These points man- 
ifest as spikes in the direct neighbor distance function, 
£(£). After a spike, there is no guarantee that the path 
is still physical. Therefore, we segment the time-series 
based on the spikes. The second view (from the second 
camera) treated with the same algorithm contains the 
information to correct such wrong segments. From the 
camera calibration the rotation matrix which transfroms 
the orientations seen by one camera into the CS of the 
other one, is known. Therefore, both views are expressed 
in an intermediate, common CS where the segments 
with d(# caml , ^ cam2 ) > 30° can be corrected. 

The algorithm presented so far assumes an ortho- 
graphic view. This condition holds only true if the par- 
ticle center is on the optical axis of the camera or in the 
case one uses tele-centric lenses. In the present exper- 
iment we do not, and the perspective effect alters the 
measured orientation (note that the parallax displace- 
ment corresponds to a change in the 2D projection, and 
hence to a rotation). The distortion induced by the per- 
spective is characterized by the position of the particle 
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center in the camera image, X, and the focal length, 
/. Common camera objectives allow only small angles, 
Tpersp = atan (||X|| //) < 15°. As a consequence we 
assume that the shape of the particle does not change 
and we introduce an orientation matrix R pe rsp. (taking 
advantage of the Rodrigues formula Eq. (9)): 



R 



(-y,x,o) 



, at an 




such that the measured orientation is related to the 
absolute orientation ahs by R = IL persp £L(# a bs)- The 
perspective distortion can then be removed from the 
orientation time series. 

Finally, after correcting for perspective distortion, a 
combined time-series of orientation can be built using the 
information from both views, if they are expressed in the 
same CS. Euler angles are not locally orthogonal, hence, 
we use the weighted mean of the orientation expressed 
in the axis-angle representation. The variance within a 
moving window of the direct neighbor distance function, 
£(£), proves to be a good error estimator of the noise, 
since for short times the particle is assumed to rotate 
smoothly. A sample orientation track is shown in the 
upper panel of Fig. 7. 



360 




0.4 0.6 
time [s] 



FIG. 7: A sample orientation track; it is X = Q, y = +, 
9 Z = □, the bottom plot shows the distance (in degrees) be- 
tween the independent orientation measurements from the 2 
cameras. 



particle texture must be first tuned in order to optimize 
these parameters - by trial and error methods. For the 
orientation algorithm per se, we have used a series of syn- 
thetic images of known orientation. We found that the 
measurement error is 2° , which is smaller than the size, 
Afine = 3° of the fine grid used in the image process- 
ing (cf. paragraph IIIB lc). A finer grid would improve 
the resolution for ideal images, but not for real images 
which, as stated above, always contain some amount of 
distortions or impurities. In addition, the fast dynamics 
of the particle and high frame rate ensure that wrong 
detection do not persist for longer than a few frames. As 
a result, most defects are detected and skipped or in- 
terpolated or handled as part of post-processing (wrong 
orientations correspond to jumps in the direct neighbor 
distance function). 

We illustrate the accuracy of the detection on 2 ex- 
amples. The first one concerns the agreement between 
the orientation as estimated from each camera measure- 
ment. In the upper panel of Fig. 7, the combined 3 angles 
with respect to the Lab coordinate system are plotted. 
The lower panel shows the distance (in degrees of an- 
gle) between the two estimations, d (# cam i, CSLUl2 )- The 
probability density function (PDF) of these distances, 
computed with and without processing for perspective 
corrections are shown in Fig. 8. When the correction for 
perspective distortion is made, a weighted average leads 
to an absolute error equal to 3.5°. 



Perspective removed 
x Raw data 




5 10 15 

d(e cam1 ,e cam2 ) [deg] 

FIG. 8: Probability density function (PDF) of the distance 
between the orientations measured from cameras 1 and 2, 
without correction for perspective distortion (x) and with it 

(°)- 



C. Robustness 

A full study of the accuracy and robustness consider- 
ing all possible distortions is beyond the scope of this 
article. In practice, the problems with real images are 
mainly caused by reflections, bad illumination, and ob- 
jects (such as bubbles or dirt particles) between the par- 
ticle and the camera. The setup, light conditions and 



IV. RESULTS 

The results in this section correspond to the flow cre- 
ated by counter-rotation of the driving disks at a rate of 
3Hz. In this case the power injection is of the order of 
e ~ 1.7 W/kg, the integral time scale Tj, is about 0.3 s, so 
that the dissipative time and space scales are r] ~ 30 fim 
and t v ~ 1 ms. As a result, the particle tracked has a 
size corresponding to D/rj ~ 600 and D/L[ nt ~ 0.6 (Lint 
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FIG. 9: Histogram of track segments. The exponential decay 
rate is of the order of the integral time Tl of the flow. 
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is the scale at which energy is fed into the flow). The 
flow Reynolds number based on the Taylor micro-scale 
is R\ ~ 300. The camera frame-rate is 600 Hz, and the 
trajectories analyzed have been selected so that their du- 
ration is longer than 0.25 Tl and most range between 0.5 
and 3 Tl. 



TABLE I: Characteristic values (mean ± rms) for the parti- 
cle motion. The angular variables are given for the lab and 
particle coordinate systems. 





X 


y 


z 


Norm 


V 


[m/s] 


0±0.28 


± 0.40 


0±0.37 


0.6 ±0.2 


a 


[m/s 2 ] 


-0.1 ±5.5 


-0.3 ±5.6 


0.2 ±6.0 


8.4 ±5.3 


w L 
w p 


[rad/s] 
[rad/s] 


-0.1 ±8.0 
-0.1 ±7.1 


0.2 ±7.8 
-0.3 ±8.2 


0.2 ±7.7 
-0.1 ±8.1 


12.2 ±5.8 


a L 
a p 


[rad/s 2 ] 
[rad/s 2 ] 


-3 ± 590 
0±624 


1±554 
2 ±516 


0±559 
1±534 


820 ± 530 



Fig. 9 shows a histogram of the duration of recorded 
tracks for which the 6D coordinates of the particle are 
recorded. It has an exponential tail (as it was also the 
case when using acoustic tracking [ ]). For very small 
times the histogram is biased by the fact than tracks 
shorter than 50 contiguous frames are discarded. Note 
also that long tracks are likely to correspond to trajecto- 
ries spanning the flow volume, i.e. a spatial extend over 
which the large scale (anisotropic) circulation cannot be 
ignored. 

However, one first result is that the particle explores 
uniformly the orientation space. This is seen in Fig. 10 
showing the probability distribution functions of the Eu- 
ler angles: as expected from a random distribution of 
orientations, the 6 X and Z components have a flat dis- 
tributions spanning a [— 7r, ±7r[ interval, with the inner 
angle y having a cos(9 y ) distribution over [— 7r/2, tt/2]. 

Interesting features are observed for the rotation dy- 
namics. The statistics of angular velocity fluctuations are 
shown in Fig. 11. The distributions are symmetric. The 
three components, with respect to the Lab CS, follow the 
same statistics. This reflects the spherical symmetry of 



FIG. 10: PDF of the orientation = {6 x ,0 y ,0 z }, the solid 
lines correspond to a uniform sampling of the orientation 
space 

the particle; furthermore, it also shows that the turbulent 
swirls at the scale of the particle have no preferred ori- 
entation. The mean of the angular velocity components 
(with respect to the Lab reference frame) is essentially 
zero, up to statistical error. The rms amplitude of angu- 
lar velocity fluctuation is of the order of 12 rad/s which 
is of the order of u rms /D = 30 rad/s. That is, it corre- 
sponds to the rotation that would result from imposing 
a velocity difference equal to almost u rms across the di- 
ameter D of the sphere. Note that it is also of the order 
of the rotation rate of the driving disks. The PDF them- 
selves displays weakly stretched-exponential tails; for a 
quantitative estimation we use the fitting function: 

e 3a 2 /2 / /ln|x/ v / 3|±2a 2 \\ , x 

"• (,) - 43-( 1 - t ( aji J J <17> 

which has been used extensively in the analysis of the 
intermittency of the translational motion of Lagrangian 
tracers [ ] - it stems from the approximation that the 
norm of the vector has a lognormal distribution. For the 
angular velocity, one finds a fitting parameter a = 0.45, 
which corresponds to a flatness factor F = 4. It would be 
F = 3 for Gaussian statistics, so that our measurements 
show only a slightly non- Gaussian behavior for the 
angular velocity. This differs from the translational 
velocity, which is found to be slightly sub-Gaussian. 

The angular acceleration has a strong non- Gaussian 
behavior, as seen in Fig. 12. Again, the three components 
follow identical statistics: there is no preferred direction 
for the torques acting on the moving sphere (with respect 
to the Lab reference frame only - the issue of lift forces 
is addressed elsewhere [ ]). The rms amplitude of an- 
gular acceleration is about 800 rad/s 2 , again of the order 
of (u rms /D) 2 . The statistics is strongly non Gaussian, 
a fit using the same stretched exponential distribution 
yields a = 0.6, i.e. a flatness factor F ~ 7.6. The an- 
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(oo-<oo >)/std(oo) [dim less] 

FIG. 11: PDF of the (normalizd) components of the angular 
velocity, uj x = O) w y — +5 00 z = □, the dotted curve is a 
Gaussian and the dashed one shows a stretched exponential 
with a = 0.45 (F = 4). 

gular acceleration can be viewed as an angular velocity 
increment over a very short time lag. Hence, the PDFs of 
angular velocity increments change shape with the length 
of the time lag - from the one in Fig. 11 for small time 
increments to the one in Fig. 12 for integral times. 

For comparison, we recall some features of the transla- 
tional dynamics of the particle. It has statistical charac- 
teristics which are very close to the one reported for neu- 
trally buoyant inertial particles with a size much closer 
to the dissipation scales of turbulence [3, 9-11]. The 
translational velocity follows a Gaussian distribution, its 
acceleration is strongly non-Gaussian, with stretched ex- 
ponential tails. Using the stretched exponential distribu- 
tion leads to a = 0.6. One thus observes that the angular 
variables have intermittent dynamics, just as the transla- 
tional motion. The fact that it is quite pronounced, even 
for an object of size close to the integral scale of motions 
came as a surprise and deserves further investigations. 

V. CONCLUDING REMARKS 

The focus of the work reported here has been to es- 
tablish a technique for the study of angular and transla- 
tional motion of a particle freely advected by a turbulent 
flow. We have shown that the measurement technique 
is robust, efficient, and accurate. As an application, we 
report here the first observation of intermittency for the 
rotational dynamics of an inertial particle. 

We note that the algorithm used to compute the an- 
gular velocity can be applied to a set of particle attached 
to a rigid body which are tracked using standard parti- 
cle tracking algorithms. If one records the positions in 
space of 3 or more points, Pi . . . P/v at time t and t + At, 
their motion can be split up into a translation of their 



center of mass (CM) plus a rotation. Once the trans- 
lation part is subtracted, the rotation, £ kabsch ? of the 
points Pi . . . P/v around their CM can be computed effi- 




-10 -5 5 10 



(a- <a» / std(a) [dim less] 

FIG. 12: PDF of angular acceleration; it is ol x = cx y = +, 
ol z = □, the dotted curve is a Gaussian and the dashed one 
shows a stretched exponential with a = 0.6 (F — 7.6) 



ciently using Kabsch's [26, 27] algorithm. S kabsch is then 
the matrix representation of the change in orientation, 
and the angular velocity, u> F , (in the particle reference 
frame) at time t can be extracted as done here. It should 
be pointed out that, one does not gain access to neither 
the angular velocity in the Lab reference frame, u; L , nor 
to the absolute orientation, 6. 

The strong intermittency in the particle's rotation may 
eventually be traced back to the complex interaction be- 
tween the particle and its wake. One notes that this 
is inherently a finite size effect; for particles with very 
small diameters (compared to the Kolmogorov length) 
the translational and rotational dynamics are note cou- 
pled. For larger particles, as in our case, the influence 
of rotation on the motion of the particle is of interest, 
and will be the object of further analysis. One may 
also note that the influence of the inhomogeneity at large 
scale must be clarified. Further measurements in a more 
isotropic turbulent flow (such as the Lagrangian Explo- 
ration Module [ ]) are underway. 
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